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ABSTRACT 

A new formalism is presented for finding equilibrium distribution functions for axisym- 
metric systems. The formalism, obtainded by using the conce pt o f fractional deriva- 
tives, generalizes the methods of lFrickd (|1952| ). iKalnaisI (|1972| ) and lJiang fc Ossipkovl 
(2007;), and has the advantage that can be applied to a wider variety of models. We 
found that this approach can be applied both to tridimensional systems and to flat 
systems, without the necessity of dealing with a pseudo- volume mass density. As an 
application, we obtain the distribution functions of the Binney's logarithmic model 
and of the Mestel disc. 
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1 INTRODUCTION 

The construction of self-consistent models for stellar systems 
is of great interest in astrophysics. As it was pointed out 
by some authors, the most straightforward way to perform 
such models is to start with an assumed potential defining 
the mass density p (via Poisson's equation) and the families 
of orbits that can lie within the system (via Newton's equa- 
tions of motion). Since p is the integration of the distribution 
function (DF) over the velocity variables in the phase space 
of the system, the problem of finding the DF is that of solv- 
ing an integral equation. For that reason, such procedure 
is called as th e "from p to /" approach f or finding a self- 
consis t ent DF | Binncy & Trcmainel l| 19871 ) : lHunter fc QuianI 
(|l993h : iJiang fc Ossipkovl (|2007l )) and the integral relation 
connecting / and p is known as the self-gravitation equation. 

By Jeans's theorem, an equilibrium DF is a function of 
the isolating integrals of motion that are conserved in each 
orbit. It has been shown that, for certain potential-density 
pairs (PDP), it is possible to find analytically such kind of 
DFs. The simplest case of physical interest corresponds to 
spherically symmetric PDP, described by an isotropic DF 
that only depends on the energy. Eddington (1916) showed 
that it is possible to obtain such DFs by first expressing the 
density as a function of the potential, and then solving an 
Abel integral equation. 

Another case of interest in astrophysics is the corre- 
sponding to axially symmetric systems, for which a great 
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variety o f PDP has been construc t ed, e.g. IKuzmin (1 19561): 
iToomrel il 19631): liyrivamotcl dl97lh: iBaginl 7l987l); iKalnais 
J 19721) : iMivamoto fc Nagail (|l975T ); iMivamoto fc Nagai 
jl976l i: Kutuzov & Ossipkov (1980, 1986, 1988); Evans 
(1993, 1994); Kutuzov (1995); J iang (2000); Jiang & Moss 
(2002): lGonzalez fc Reinal (|2006l ): Ossipkov & Jiang (2007). 
Now then, the first approa c h to obtain axisymmetric DFs 
was introduced by iFrickel (|l952l ) and, after that, differ- 
ent integral transformation techniques were used to obtain 
two-integral DFs (e.g. Lynden-Bell (1962); Hunter (1975); 
Kalnajs (1976); Dejonghe (1986); Hunter & Qian (1993)) 
but there appear several obstacles concerning to the validity 
of such transformations and the requirement of the complex 
analyticity of a density-related integral kernel to complex 
arguments. 

Recently, Ijiang fc Ossipkovl (|2007l ) presented a new 
method for the axially symmetric case, where the equilib- 
rium DF depends on the energy E and the angular mo- 
mentum about the axis of symmetry L z , i.e. the two clas- 
sical integrals of motion. They developed a formalism that 
essentially combines both the Eddington formulae and the 
Fricke (1952) expansion in order to obtain the DF's even 
part, starting from a density that can be expressed as a 
function of the radial coordinate and the gravitational po- 
tential. Thus, for a given p, the corresponding even DF can 
be obtained by solving an Abel integral equation. Once such 
even part is determined, one can find the DF's odd part by 
introducing some reasonable assumptions about the mean 
circular velocity or using the maximum entropy principle. 

In the present pap er we show a n ex tension of the for- 
mulae derived earlier bv lFrickel (|l952l ) and lJiang fc Ossipkovl 
l|2007l ). by introducing the fractional derivative concept. The 
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approach developed here has several advantages over the 
methods introduced before. In one hand, the mathematical 
difficulties involved in the formalisms based on transforma- 
tion techniques can be easily overcome. I t is worth to point 
out t hat the formalism introduced by Ijiang fc Ossipkovl 
l|2007f ) demands the definition a volumetric pseudo-density, 
in order to be applicable to flat systems. On the other hand, 
our method can be applied directly to the case of flat sys- 
tems. 

Assume that <E> and E are, respectively, the gravitational 
potential and the energy of a star in a stellar system. It is 
useful to define a relative potential ^ = — $ + $o and a rela- 
tive energy e = — £ + $0, in su ch a way that the system h as 
only stars with energy s > jBinnev fc Tremainel (119871 )). 
For the case of an axially symmetric system, it is custom- 
ary to use cylindrical polar coordinates (R, ip, z), where the 
velocity is denoted by v = (vr, v v , v z ). Such system ad- 
mits two isolating integrals: the component of the angular 
momentum about the z-axis, L z = Rv v , and the relative en- 
ergy e. Hence, the DF of a steady-state stellar system in an 
axisymmetric potential can be expressed as a non-negative 
function of e and L z , that vanishes for e < 0, denoted by 
f(e,L z ) and related to the mass density as 



4yr J-V j-Ry/HV-e) 



R 



f+(e,L z )dL z de, (1) 



where /+ (e, L z ) is the even part of / with respect to the an- 
gular momentum L z . On the other hand, for the case of flat 
systems, the surface mass density S is related to / through 



s = 4 /"* f f+ (e,L z )dL z ds 

Jo Jo ^2R 2 (V-e)-L 2 ' 

Now, by defining a pseudo-volume density p, according to 
(Hunter and Quian (1993)) 



■ = V2 



* S(i? 2 ,^ / )dv&' 



(3) 



it is also possible to use equation (JTJ to deal with these flat 
systems. 



where ct n S R and the 2 in the exponent of L z guarantees 
that /+ is even. Now, performing the integral {1} with re- 
spect to L z , we obtain 



E 
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Me)(* - 



' 2 de, 



(5) 



for a n > —1/2, while it diverges for a n < —1/2. Therefore, 
we assume that the the mass density is given by 



p( J R,*) = ^ J R 2 ""p„(*), for a n >-\, 



(6) 



A comparison between ([5} and ([6]), leads to the relation 



Pn(*) 
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(7) 



At this point, we introduce the fractional derivative operator 
D 3 X , which represents a 7-ord er derivative, with re s pect to 
x, for any real value of j (see lBologna fc Grigolinil l|2003t) ). 
Assuming that (£>^/5n(*))*=o = for all j £ (0, a„ + 1/2), 
then 



h n {e)de. (8) 



This integral equation is simpler than the first one and can 
be inverted easily if one takes the derivative once again with 
respect to ^, 



h n {e) 



7r2 Q "+2r(a„ + i) 
and the distribution function can be expressed as 



f+(e,L z )=J2- 



5 Pn(*) 
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tr(a n + |) 



(9) 



(10) 



When a n € N, by the definition of the Riemann-Liouville 
operator, equation fllOD reduces to the formulae obtained by 
Ijiang fc Ossipkovl \200H) . 

As a particular case, suppose that fini^) can be written 
in the form 



0k 



(11) 



2 TRIDIMENSIONAL SYSTEMS 

Most of the methods developed to solve require some 
kind of dependence between the DF, t he rela t ive e nergy 
g and the angular mom entum L z (see iFrickd (jl952l ) and 
Ijiang fc Ossipkovl l|2007h . as examples) , which will define the 
corresponding mass density p(_R, \&). Therefore, in order to 
study the problem, we shall start by assuming some general 
types of DFs. 



2.1 DFs of the form £ L 2 z a " h n (e 



So, applying the fractional derivative operator to (|ll|l 

Y[p h -a n - 2) 
for ft > o„ + 1/2, and the corresponding DF is 

A nk T(l3 k + l)L 2 z a -e^- a "-i 



t 7r2""+2r(a„ + |)r(ft - o„ - \) 



(12) 



(13) 



This result is totally equivalent to the Fricke solution, for 
real values of a n , and therefore can be considered as a gen- 
eralization. 



To generalize the Jiang & Ossipkov method, we suppose that 
the DF can be written as 



f + (e,L z ) = Y t Ll an hn(e), 



(4) 



2.2 DFs of the form ^ L 2 z arl g n {Q) 

n 

It is possible to derive a more general expression for the DF 
if we assume that it depends on e through Q = e~L 2 /(2R 2 ), 
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where R a is a scaling radius. Suppose that the system has 
only stars with Q > 0, so f(Q, L z ) = for Q < 0. Here, Q — > 
e as 7?a — * oo. The fundamental equation can be written, in 
terms of Q, as 



and so we obtain 



.Sl) 

l'= — I I 1 f+(Q,L z )dL z dQ, (14) 

'o io 

where f+(Q,L z ) is the even part of f(Q,L z ). So, following 
a similar procedure than in section 12.11 one can find that a 
DF of the form 



f+(Q,L z ) = J2- 



D^ + *p n (V) 



* = Q 



7r2""+2r(a n + i; 

corresponds to a mass density of the form 

i? 2Q "Pn(*) 



(1 + - 

n \ ± T H : 



-R 2 W + J 



(15) 



(16) 



for q„ > —1/2, where a n G R. Now, if we sum over all 
posible values of R a we obtain the general expression 



f+(Q,L z ) = J2' 



* = Q 



(17) 



(18) 



„,. 7r2""+2r(a„ + i) 
corresponding to a density of the form 

with i? a > and a n > —1/2. 



2.3 Models with divergent gravitational potential 

In a system in which the gravitational potential has no up- 
per bound, it is not possible to define correctly the relative 
potential ^ and the relative energy e, because the scape en- 
ergy of the system is oo. For this reason, we shall write the 
fundamental equation in terms of E and $, 



p(R,*) = 



f+(E,L z )dL z dE, (19) 
an be wi 

f+(E,L z )=J2LT"h n (E), (20) 



/* Jo 

and we will suppose that the DF can be written as 



for a n > —1/2, and that the density is given by 

(21) 

n 

for a n > —1/2. Thus then, by integrating with respect to 
L,, follows that 



4tt2 



2a n + ± JiS> 
Now, if we assume that 



(E)(E -§) an + zdE. (22) 



lim D|p„($) = 



(23) 



for all j € (0, a n + 1/2), then 



i -I poo 

Dl n + *p n =J7r(-2) a "25r(a n + -)y h n (E)dE, (24) 



h n (E) 









$ = E 



7r(-2)-"+§r(a„ + i; 



Therefore, the distribution function is 



D 



? Pn($) 



f+(E, L z ) = V 3 

V 7r(-2)-+§r(an + |) 



(25) 



(26) 



for a n > -1/2. 

On the other hand, we can assume that the density has 
the more general form 

i? 2 ""Pn(<&) 



for a n > —1/2. So, the corresponding DF will be 



(27) 



L 2 , 



f+(Q,L z )=Y J - 



? Pn($) 



*=Q 



(28) 



Jr (_2)-»+*r(o n + i) 

for q„ > -1/2 and Q defined as Q = £ + L 2 /(2i? 2 ). Now, 
if we sum over all the posible values of R a , we can obtain 
the generalization 



f+(Q,L z ) = J2' 



l> Dir*~ Pn {$>) 



~ 7r(-2)"" + 5r(a n + I) 
corresponding to a density of the form 

i? 2Q "Pn($) 



p(7?,$) = ^ 



with i? a > and a n > -1/2. 



(1 + 



(29) 



(30) 



3 FLAT SYSTEMS 

The formalism sketched above can also be used directly in 
th e case of flat system s . Not e that in the method introduced 
by IJiang fc Ossipkovl (|2007T ) it was not possible, since the 
fundamental equation could not be solved using the Abel 
integral equation. Now, we will proceed similarly to the 
tridimensional case, finding the DFs for different kinds of 
densities. Then, we will also study the case of models with 
divergent gravitational potential. 



3.1 DFs of the form ^ L 2 z c " n h n {e) 

n 

As in the tridimensional case, first we suppose that 

f + (e,L z )=J2 L2 " n M £ )' (31) 

n 

So, by integrating (J2j) with respect to L 2 , we obtain 

S - E f<TZt n + \l A* - (32) 



2-(«n+Dr(a„ + 1) J 
for a n > —1/2. Therefore, if we assume that 

E(fl,*) = ^i? 2a " ( T n (*) ) (33) 
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a„ > —1/2, it follows that 
V^T(a n + §) 



<t„($) = 



2-(««+ 1 )r(a„ + 1) y 



(* -e) an h n {e)de. (34) 



Now, if (-D^crn(v&))*=o = for all j € (0, a„), then equation 
(1341) leads to 



(35) 



D*"<t„(*) = ^2 a " + 1 r(a„ + I) /* h n (e)de. 

2 Jo 

Consequently, 

£>S n+1 er n (#)| T 
kn{£) - 0F2-+iI> n + I)' 



and the DF is 

f+(e,L z ) = J2 



0F2-"+ 1 r(a r , 



(36) 



(37) 



This equa tion, for a n = 0, corresponds to the method de- 
veloped in Kalnais (1976), working in an adequate rotating 
frame. 

As a particular case, suppose that 



Then, by taking the fractional derivative we obtain 

^" + V4*) = E- fcr(/3fc + 1) *' 5fc " Q "" 1 

fc 

and the DF is 



f+(e,L*)=J2 



r(/3 fe - a n ) 

0F2-+ 1 r(« n + ±)r(/fc - «„) 



(38) 



(39) 



(40) 



This relation can be interpreted as the analogous case of the 
Fricke expansion, when we are dealing with flat systems. It 
can be verified performing the pseudo-volume density © of 
R 2a " ^ /3k and taking the Frike component corresponding to 
the tridimensional case. 



3.2 Other DFs for flat systems 

We can generalize the result (|37jl if we express the DF in 
terms of Q — e — L 2 /(2R 2 ). In this way, if the density has 
the form 



iiu + fr)"- 



the corresponding DF is 



f+(Q,L z ) = J2 



0F2«»+ 1 r(a„ + §) 



(41) 



(42) 



for R a > and a„ > -1/2. 

Furthermore, if we consider models with gravitational 
potential having no upper bound, we can deduce that for a 
density 



Z. (1 + #)an+ i 



(43) 



and assuming that lim^^^ _D| ) cr„($) = for all j G (0, ct n ), 
then 



L 2„ n ££»+l ffn ($J 



?r, (44) 



( n V5F(-2)C"+ 1 )r(a„ + 2 , 
for i? a > 0, a n > -1/2 and Q = £ + L 2 z /(2R 2 a ). 

4 SOME APPLICATIONS 

In this section we will use the formulae introduced above 
to the Binney's logarithmic model and the Mestel disc, and 
we will see that their corresponding DFs match exactly with 
those that were found through the application of other meth- 
ods. Binney's logarithmic model has a gravitational poten- 
tial of the form 



$(R,z) = h 2 \n(l + R 2 + ^ ) Mo) 



whereas its mass density is 



p(R,z) 



^(l + 2g 2 + J R 2 + (2-g- 



47rG<? 2 (l + R 2 + z 2 q~ 



that can be written as 



^ 2 {[(l-g 2 )fl 2 + l] + ( g 2 -i) e 2 */"o 2 } 



(46) 



(47) 



27rGg 2 e 4 */"o 
So, as D"e ax = a a e ax for any a € R, we obtain 

f+(E,L z ) = AL 2 z e' 4E/v ° + Be~ iE/v2 o + Ce' 2E/v °, (48) 
where 



y/2 



G = 



2q 2 



7rfGq 2 u§ ' " TT§Gq 2 v ' ^l 2 Gq 2 VQ 

the same DF founded in Ijiang fc Ossipkovl J2007h . using 
the Abel's integral equation, and in lEvansl l|l993bT ) using 
Lynden-Bell's method. 

Another case of interest is the Mestel disc, characterized 
by a gravitational potential of the form 

R 



§{R) = v 2 c In 



Ra 



and a surface mass density given by 

Ro 



B(Jl) = Eo- 



R 



(49) 



(50) 



where Eo = v 2 /(2ivGRo)- Now, we can write (|50|) as 

E( J R) = jjam Eo.e-Cam+D*/.^ (51) 

i? m 



for any m € R. So, equation (|44|) for R a — > oo leads to 

/+(£,£,) = FL 2m e- E/CT2 , (52) 
where F and a are the constants given by 

E 2-( m + 1 ^- 1 / 2 



-, F: 



(53) 



2m + 1' " r(m + l/2)i?g m cr 2m + 2 ' 

This solution was obtai ned as well by lEvansI |l993b) , and 
was proposed earlier bv lToomrel (|l977h . 
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5 DISCUSSION Toomre, A., 1963, Ap. J., 138, 385. 

. , , , , , , , Toomre, A., 1977, ARA&A, 15, 437. 

In contrast with the methods based on integral transfor- 
mation techniques, the formalism developed here does not 
require that the mass density has an analytic continuation 
to complex arguments. Indeed, this is the principal disad- 
vantage involved in such methods. Moreover, our fractional 
derivative approach can be regarded as a general method 
that contain s, as part i cular case s, the results obtained b y 
iFrickd i|l952D , lKalnaisl (|l976l ) and I Jiang fc Ossipkovl (|2007h . 
The method developed here can be applied to a wider va- 
riety of axisymmetric models, due to the generic form of 
the density as a function. Another advantage of this formal- 
ism is that it can be applied directly both to tridimensional 
systems and to flat systems, without the implementation 
of a pseudo-volume density. Therefore, taking into account 
all the above statements, the present formalism represents 
a powerful tool on the construction of self-consistent stellar 
models. 
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